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We present a new methodology to solve the Anderson impurity model, in the context of dynamical 
mean-field theory, based on the exact diagonalization method. We propose a strategy to effectively 
refine the exact diagonalization solver by combining a finite-temperature Lanczos algorithm with an 
adapted version of the cluster perturbation theory. We show that the augmented diagonalization 
yields an improved accuracy in the description of the spectral function of the single-band Hubbard 
model and is a reliable approach for a full d-orbital manifold calculation. 



The understanding of materials with strongly corre- 
lated electrons is one of the main challenges of mod- 
ern solid state physics. Triggered by the discovery 
of high-temperature superconductivity in copper-oxides, 
the study of doped Mott insulators has grown in the last 
decades, building on the development of theoretical tools 
designed to solve accurately models of strongly correlated 
electrons. 

Despite that the exact solution of simple correlated 
theoretical models in two or three dimensions is lacking, 
accurate predictions for the properties of strongly cor- 
related solids are obtained by using approximations [1]. 
A central role has been played by the dynamical mean 
field theory (DMFT)[2], a nonperturbative method which 
allowed for the first complete description of the Mott- 
Hubbard transition. This method has been extended to 
a variety of correlated methods and combined with den- 
sity functional theory, leading to remarkable agreement 
with the properties of many correlated materials. 

Within DMFT the original lattice problem is mapped 
onto an effective Anderson Impurity Model (AIM) , which 
can be solved numerically. The numerical solution of the 
AIM is the real bottleneck of a DMFT calculations, and 
the development of impurity solvers which are both ac- 
curate and computational-wise efficient is a very active 
line of research. Among the many solvers proposed and 
extensively tested through the last decade, those based 
on exact diagonalization (ED) [3-5] or on the continuous- 
time quantum Monte Carlo (CTQMC) [6] are used ex- 
tensively in this line of research. 

Indeed, the recent development of CTQMC has gener- 
ated a strong activity in the field. For single-site DMFT, 
CTQMC yields an exact solution of the AIM problem 
within the statistical error bars in imaginary time. The 
main limitation of the approach is that the evaluation 
of real-frequency spectra requires a poorly conditioned 
analytical continuation, based on the maximum entropy 
method [7] (or some alternative strategy), which strongly 
limits the possibility to study fine details of the spec- 
tra. For multi-orbital or cluster extensions of DMFT, 
CTQMC suffers from the fermionic sign-problem as long 
as inter-orbital hybridizations are present, and it is there- 



fore limited to finite temperatures. 

The ED solvers arc instead based on a finite discretiza- 
tion of the AIM, through the representation of the effec- 
tive bath in terms of a small number of bath-sites. In 
practical implementations the bath size (A^) is severely 
limited because of the growth of the Hilbert space: its 
dimension scales exponentially with the total number of 
sites N s (bath sites and impurity orbitals). Nonethe- 
less, the use of Lanczos-based algorithms allows to deal 
with large Hilbert spaces, and the discretization at low 
temperature[4, 8] is fine enough to accurately compute 
the thermodynamic and static observables. In particu- 
lar, the finite size effects affect the spectral functions, 
which arc slowly converging to the continuous features of 
the exact DMFT solution. Notwithstanding reasonably 
accurate static phase diagrams obtained with the Lanc- 
zos solver for three[5, 9] and five[10] orbitals system, the 
limitation in the bath size becomes also particularly rel- 
evant for multi-orbital AIM models, which are necessary 
to realistically describe transition metal oxides, as the 
impurity's orbitals (five for a d manyfold) contribute to 
the enlargement of the Hilbert space. 

In the view of the increasing interest in the effect 
of multi-orbital and multi-site Coulomb correlations in 
transition metal oxides, an improvement of the ED 
methodology is hence called for. 

In this work we propose a simple, yet efficient, mod- 
ification of the ED solver for the DMFT, which allows 
to treat larger bath sizes without increasing the compu- 
tational cost. In particular, we propose a combination 
of ED, which is used to solve a small part of the bath, 
but we include an extra set of bath levels which are per- 
turbatively coupled to the first part using concepts along 
the line of cluster perturbation theory (CPT) [11, 12]. 
The augmented effective bath in the ED-CPT method 
improves the accuracy in both static and spectral prop- 
erties as compared with standard ED, and it allows for 
accurate calculations for multi-orbital models. 

Within DMFT, the effective AIM is subject to a self- 
consistency condition which relates the Green's function 
of the impurity model G(iuj n ) to the so-called Weiss field 
Q^^iojn), which completely characterizes the AIM. For 
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FIG. 1. (Color online) Sketch of the AIM in ED-CPT: the 
impurity (center circle) is connected to JV& nearest neighbor 
bath sites (stars), and each of the latter is connected by a 
small perturbative parameter Vcpt to a line of additional 
sites (squares, N c sites in each chain). The inner system (light 
filling) is solved by using a standard Lanczos algorithm. The 
Green's function of the full system (inner and outer region) is 
obtained by treating the full system by cluster perturbation 
theory. 

the single band case, the Hamiltonian of the AIM reads: 
H = ^ e p<jal a a p(J +^2 Vpa{a\, a d a + hc)+Un^ni (1) 
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where at [a pa ) creates (destroys) a particle with spin 
a in the p-orbitals (p € [1, N p ]) of the bath and d) a (da) 
creates (destroys) a spin a particle on the impurity, and 
U is the static Coulomb repulsion on the impurity. Spin 
indices are omitted throughout this letter. For a fixed 
set of the Hamiltonian parameters, we solve the AIM (1) 
by using a Lanczos algorithm to converge the low-lying 
states of the spectrum[4, 5] which contribute to the ther- 
mal average: we consider an energy cutoff E max such that 
the Boltzman weight of the discarded states is negligible 
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< 0.001, where E is the ground state en- 
ergy. Once the eigenstates are obtained we compute the 
dynamical and static observables[13]. 

In the DMFT methodology, the Weiss field Q^ 1 (iw n ) = 
iuj n + fj, — A(iu) n ), obtained through the mean- field self- 
consistency condition, must be represented in the discrete 
basis of the AIM. To this intent, we minimize a suitable 
distance between the Weiss field (or, equivalently, the 
function A) and the discretized AIM hybridization func- 



tion A ED {iuj n ) = YfU V*/{iw n - e 



d = 
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\A ED (ioj n ) - A (iu> n )\ K, 



(2) 



where ujq is a hard cutoff on the number of Matsubara 
frequencies (we set uj = 100 throughout the rest of this 
work). The denominator is a weighting factor that can be 
used to improve the precision on the low energy frequen- 
cies close to the metal-insulator transition [5, 14]. The 
distance d is furthermore renormalized by the number of 
frequencies uj < loq. 

The minimization (2) boils down to the problem of 
fitting an arbitrarily complicated function with a finite 
number of rational functions. This is the only source of 
potential systematic errors and limitations of the stan- 
dard ED algorithm. Indeed, for low temperatures the 
number of bath-orbitals required to resolve the low en- 
ergy features characterizing the hybridization function 
can be large. Noteworthy, there is no simple limit where 
the fitting method used in the ED solver is exact, in 
particular the finite size effects and systematic errors in 
the fitting procedure are equally severe for the trivial 
U = and large U limits. Moreover, the fitting pro- 
cedure may become problematic near the zero repulsion 
U ~ 0, where the hybridization might have a complicated 
structure difficult to capture with the fitting method. It 
is worth reminding that, while there is no strong limi- 
tation to the number of bath-orbitals at the level of the 
fit, the size of the Hilbert space of the effective AIM is 
exponentially growing with the number of sites, rapidly 
becoming out of reach for any ED-based algorithm. In- 
deed, although the Lanczos solver yields static and dy- 
namical observables with a remarkable precision, the ef- 
fective discretized AIM might very poorly reproduce the 
DMFT hybridization, giving in turn large errors in the 
spectral functions. 

In this work, we solve this problem by extending the 
AIM to a larger system. This extension is realized by 
coupling additional lines of N c non-interacting sites to 
each of the Nb bath-orbitals (see Fig. 1). In the spirit 
of the CPT, the couplings V p ^ C pt between each chain and 
the corresponding bath-site are introduced as small per- 
turbative parameters[12]. Indeed, we treat the hybridiza- 
tion between the AIM and the chains within CPT, which 
amounts to use the following expression for the impurity 
Green's function: 



iCPT 



- V 



CPT 



(3) 



where G CPT (iui n ) is the full Green's function of the cou- 
pled (AIM+chains) system, G Vo (iu) n ) is a block-diagonal 
matrix which contains the Green's function of the decou- 
pled AIM and chains, and V CPT is the hybridization ma- 
trix that connects the two subsystems. All the matrices 
are defined in the space composed by the bath-orbitals 
and the chains. The hybridization function of the ED- 
CPT explicitly depends on the additional CPT parame- 
ters, and can be written in orbital space as A ED (iu> n ) 
V (iui n — e)~ V [15] , where V and e are respectively the 
matrices of the impurity-bath and bath-bath couplings. 
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FIG. 2. (Color online) Spectral function p(u>) obtained for 
the Hubbard model solved by ED-CPT at a) U/D = 1.5 and 
c) U/D — 4. For comparison, the spectral functions obtained 
by ED at b) U/D = 1.5 and d) U/D = 4 are shown. The ED 
results with iV;, = 3 (diamonds) and N b = 7 (squares) bath 
sites are compared with the DMRG results (circles). The 
CPT-ED results obtained with N b = 3 inner bath sites and 
augmented with lines of N c = 6 sites (diamonds) and respec- 
tively Nt = 7 and N c — 4 (squares) compare remarkably well 
to the DMRG data (circles). 

Thus, using this simple strategy, we are able to effectively 
increase the size of the bath and to improve our fitting 
procedure without increasing the size of the Hilbert space 
in the Lanczos diagonalization. The CPT is exact both 
in the [7 = limit and in the infinite-?/ limit (in which 
all the hybridizations are negligible) [12], while for inter- 
mediate value of the interaction, the method is expected 
to be reliable provided that Vqpt remains small. Thus, 
we add a simple Lagrange term A (j2 p |V^ CPT | 2 ^ to the 
fitting distance d and impose an upper hard-cutoff on 
the coupling V^ CPT [16]. We emphasize that the present 
method is easily implemented, since it only differs from 
the standard ED algorithm in the fitting procedure. 

To begin with, we benchmark the ED-CPT method 
for the single-band half-filled Hubbard model against the 
standard ED algorithm and a density matrix renormal- 
ization group (DMRG) calculations [17]. For sake of sim- 
plicity, we work on the Bethe lattice with a semi-elliptical 
density of states p(e) = ttD 2 \/D 2 — e 2 /2. For this sys- 
tem the self-consistency condition simplifies to a linear 
relation A(iu n ) = G(iu n ) (D 2 /A) [18]. 

We first compare the spectral function p(u>) = 
-ImG(w)/vT obtained by ED-CPT (see Fig. 2.a,c) with 
the one obtained using ED with the same number of 
bath orbitals (see Fig. 2.b,d). The results are bench- 
marked against DMRG calculations performed using a 
large number of bath sites N b = 30 (U = 1.5) or 
N b = 60 (U = 4). As expected, for a small num- 
ber of bath sites Nf, = 3 we observe strong deviations 
for |w|/-D < 0.2 in the ED results at small repulsion 
U/D = 1.5 (Fig. 2. b, diamonds). For this interaction, 




FIG. 3. (Color online) a) Fitting distance d obtained by 
using the ED and ED-CPT. b) Integrated absolute value of 
the difference in the spectral functions for the Hubbard model 
with U/D = 1.5. 

the system is in a metallic state, and hence a dense dis- 
cretization of the hybridization is required to recover the 
right low-energy spectral properties. This is dramati- 
cally improved in ED-CPT (Fig. 2.a,diamonds), where 
the spectral function with Nf, = 3 is essentially identical 
to the DMRG results in the low-energy range. Larger- 
energy features, such as the peaks at \uj\/D — ±0.35 are 
also improved by ED-CPT. For a larger number of bath- 
orbitals, N b = 7, both ED and ED-CPT give results close 
to the DMRG. Noteworthy, the Lagrange parameter A 
suppresses the coupling with the chain Vcpt when the 
number of bath orbitals is already sufficient to obtain a 
faithful representation of the hybridization function. Al- 
though the ED-CPT method is expected to provide a 
significant improvement for small U, we now explore the 
strong interaction regime U / D = 4 which corresponds to 
a Mott insulating solution (see Fig. 2.c,d). The finite 
discretization of the hybridization function leads to large 
peaks in the DOS, as seen in the ED results for N b = 3, 7 
(Fig. 2.d). ED-CPT indeed smears significantly those ar- 
tificial structures (Fig. 2.c) yielding a solution closer to 
the DMRG result. We note that also the spectral weight 
at the gap edge is quite improved by ED-CPT. 

In order to better assess the improvement of the ED- 
CPT method over the ED we show in Fig. 3. a the fitting 
distance d as a function of the effective bath size. In par- 
ticular, we find that at weak interaction U/D = 1.5 the 
distance d is an order of magnitude smaller for the ED- 
CPT than the ED method, up to discretizations N b < 8. 
This is especially interesting, since the upper limit N b = 7 
corresponds to largest system which can be carried out 
with the full ED. Thus, the ED-CPT can significantly im- 
prove also the full-diagonalization DMFT solver, allowing 
to perform reliable calculations at every temperature. 

In Fig. 3 we also report the frequency integral of the 
modulus of the difference between the spectral function 
at various values of Nf, and the result for the large value 
N b = 12, d G = f™^ \p(m,u)) - p(m = 12,u)\ du The be- 
havior of this quantity confirms the significant improve- 
ment of the ED-CPT solution with respect to the ED 
results, especially for N b < 4. Also for this quantity, for 
larger number of bath sites (N b > 8), no significant differ- 
ence between ED and ED-CPT is found for U/D = 1.5. 
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FIG. 4. (Color online) Imaginary-frequency Green's func- 
tion: Real (a,b) and imaginary part (c,d) of the Green's func- 
tion G for a multi orbital system (paramagnetic VO2). The 
data obtained by (a,c) ED-CPT (filled symbols) and (b,d) 
by ED (filled symbols) are compared with reference data ob- 
tained by CTQMC (open symbols) reproduced from Ref. [19]. 
The calculations include all five d orbitals (e 3 and t2 g ). 

The benchmark of the ED-CPT against ED and DMRG 
results for the Hubbard model demonstrates that when 
only few bath-orbitals are available to the AIM, the ED- 
CPT provides a very robust way to improve the ED- 
DMFT solution, yielding a systematic reduction of the 
finite size effects at a very moderate computational cost. 

Having benchmarked the validity of the ED-CPT al- 
gorithm against other methods in simple problems, we 
now turn our attention to multi-orbital systems. When 
applied to realistic materials calculation the finite-size ef- 
fects in the ED-DMFT solution become particularly se- 
vere. In particular only few bath sites become available 
for each impurity orbital, introducing a severe limitation 
of the method already for the case of d-orbitals. In or- 
der to test the improvements brought in by the ED-CPT 
method we apply this method to the Vanadium Diox- 
ide [20]. VO2 is a transition-metal oxide with an open 
d-shell, and hence it requires an impurity model with 
5 impurity levels [21]. This system is moderately corre- 
lated metal (Ud = 4eV) with a significant Hund's cou- 
pling J = 0.68eV. We consider a bath discretization for 
this system with Nf, = 6 (N s = 11). 

In figure Fig. 4 we compare the impurity Matsubara 
Green's function obtained by the ED-CPT (Fig. 4.a,c, 
filled symbols) and by the ED (Fig. 4.b,d, filled symbols) 
against the CTQMC results, reproduced from Ref. [19] 
(open symbols). We find with ED that the density at 
the Fermi level /o(ef) = — ImG(cj = 0)/ir of the d xz and 
d yz orbitals is quenched, and the d xy has a very large 
weight at the Fermi level (Fig. 4.d). This differs from 
the CTQMC data, which show that all the t2 9 orbitals 
are contributing to the spectral function at the Fermi 



level with a weight of the same order. The ED solver is 
hence unable to capture the relevant ingredients for this 
5 orbital system. We also tried different normalization in 
equation (2) as well as more general fitting functions [22], 
but we could not significantly improve the results. 

However, by using the extended ED-CPT Lanczos 
solver, we could recover the main features at the Fermi 
level (Fig. 4.c). In particular, despite that the ED-CPT 
slightly over-estimate the density of the d xy orbital, the 
order of the respective contribution of each orbitals is 
conserved. Another failure of the ED solver is the low 
energy frequency behavior of the real part of the Green's 
function (Fig. 4.b). Indeed, the real part of the Green's 
function vanishes for particle-hole symmetric situation, 
while its positive when the spectral weight below ep has 
largest weight and vice versa. The CTQMC data show 
that the most occupied orbital is the d xy , and the ED 
data instead occupies strongly the d xz orbital instead. 
This problem is again cured in the ED-CPT (Fig. 4. a), 
which recovers remarkably the same features as in the 
CTQMC. The advantages of ED-CPT over CTQMC is 
that the real frequency observables are readily available 
from the calculations without further approximations 
such as the analytical continuation used in CTQMC. 

In conclusion, we have presented a new solver for 
dynamical mean-field theory, which systematically im- 
proves the resolution of low energy scales. This so- 
called augmented hybrid exact- diagonalization solver is 
based on the exact-diagonalization scheme, which re- 
lies on a finite discretization of the effective bath. This 
new solver systematically reduces the finite size effects 
of the ED method by considering a system of an in- 
creased size, and by treating the additional bath-levels 
within cluster perturbation theory. While this method 
requires slight changes in the implementation of the self- 
consistency condition, it introduces no further computa- 
tional cost with respect to the ED solver. We showed 
that the ED-CPT solver provides a significant and sys- 
tematic improvement over the ED for the one-band Hub- 
bard model. We extended the calculations to a multi- 
orbital realistic material (VO2), where the number of im- 
purity orbitals needed to describe the d 5 levels prevents 
using a large number of bath sites in the ED and yields 
severe finite size effects. In this case we showed that 
the ED-CPT solver is able to qualitatively improve the 
ED results and to achieve a remarkable agreement with 
CTQMC data. The methodology described in this work 
can easily be extended to any ordered phases, such as 
magnetism or superconductivity [23]. The CTQMC code 
used in Ref. [19] was written by K. Haule[24]. C.W. was 
supported by the Swiss National Foundation for Science 
(SNFS). P.B.L is supported by the US Department of En- 
ergy under FWP 70069. M.C. and A. A. acknowledge fi- 
nancial support from the European Research Council un- 
der FP7 Starting Independent Research Grant n. 240524. 
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